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Abstract 

We present a simple discrete formula for the elastic energy of a bilayer. The formula is convenient 
for rapidly computing equilibrium configurations of actuated bilayers of general initial shapes. We 
use maps of principal curvatures and minimum-curvature direction fields to analyze configurations. 
We find good agreement between the computations and an approximate analytical solution for the 
case of a rectangular bilayer. For more general shapes (simple polyiamonds), we find a range of 
typical bending behaviors: overall bending directions along longest and shortest dimensions, inward 
bending at corners, curvature intensification near boundaries, and conical bending and partitioned 
bending zones in some cases. 
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I. INTRODUCTION 



A bilayer is a thin sheet consisting of two layers, each composed of a different material. 
When the environment external to the bilayer undergoes a change in temperature or chemical 
composition, the equilibrium strains of the two layers change by different amounts due to 
their different material properties [T] . In one common type of bilayer [2] , one layer, called the 
substrate, has an equilibrium strain of zero throughout the external environmental change. 
The other layer, called the actuated layer, has an equilibrium strain which changes from zero 
to a (possibly nonuniform) value e a . Because the layers are bonded, they undergo the same 
strain at their interface, and thus both layers cannot simultaneously be uniformly at their 
equilibrium strains. However, each layer can be brought closer to its equilibrium strain, on 
average, when the bilayer bends. Then the layer along the outer circumference is stretched 
relative to the inner layer, so the average strain in each layer is different. 

Bilayers are a widely-used technology for producing precisely-controllable motions and 
shape changes for solid bodies. A common application is the household thermostat [I], 
and newer applications include self-assembling containers [31 H] , biomedical devices [5] , and 
radio-frequency switches for wireless communications [6J. Bilayers can be considered as a 
method for making new three-dimensional structures by inducing in-plane stress. Recent 
work has used imposed stresses in deformable membranes to create wrinkled patterns [7J |8] 
and other buckled and twisted shapes (9j HO] ■ 

In this study we extend a model for thin homogeneous elastic sheets to bilayers. The 
previous model has been used to study a wide range of physical problems, including elastic 
strains due to defects in materials [TT], the crumpling of paper [12], the self-assembly of 

2 



elastic sheets under magnetic forces [13] , and the bending and buckling of spheroidal pollen 
grains under osmotic pressure [Hj. The simplicity of the model has made it very easy to 
apply to a wide range of problems. 

The present model represents the bending of a bilayer in terms of a mesh for just a single 
surface — the central surface of the substrate layer. In terms of the shape of this surface, 
we derive a discrete local formula for the strain in the central surface of the actuated layer. 
Simple formulas for the elastic energy and its gradient are then obtained, which can be 
used to find bilayer shapes as energy minima using standard optimization methods. The 
computational cost of simulating a bilayer in this way is only slightly higher than that of 
simulating a homogeneous elastic sheet, because the strain in the central surface of the 
actuated layer is proportional to the curvature of the substrate. 

The organization of the paper is as follows. In section [TTJ, we derive the bilayer model. 



In section III we present sample simulations of a rectangular bilayer, and compare with an 



approximate analytical solution. Section IV presents examples of bilayer bending for various 



bilayer shapes, most of which are simple polyiamonds. Section [V] summarizes the results. 
II. MODEL 

We represent a bilayer in discretized form using a triangular mesh. The mesh is an 
equilateral triangular lattice in the undeformed state. The substrate layer has stretching 
and bending energies which approximate those of a uniform, isotropic elastic sheet. Both 
energies may be obtained by summing simple quantities over the edges of the lattice [TT] . 
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The stretching energy is 

E s = 1 -C s J2(\r>-rj\-de q ) 2 . (1) 

Here C s is a stretching stiffness constant, d eq is the length of the edges in the undeformed 
mesh, and the sum is over distinct nearest-neighbor pairs of points. The bending energy is 

^ = a^(l-n a - n/3 ). (2) 

Cb is a bending stiffness constant, and the sum is over distinct nearest-neighbor pairs of 
triangles, with the unit normal vectors to each triangle given by n a and n^. In the un- 
deformed state, the mesh lies in the x-y plane, and all of the normal vectors are e z (i.e. 
pointing upwards in the z direction). The normal vector maintains the same direction with 
respect to its triangle during deformations. 

Seung and Nelson showed that when the stretching strain is small and the radii of curva- 
ture of the sheet are large relative to d eq , E s and Ej, converge to the stretching and bending 
energies of a uniform isotropic elastic sheet with Poisson ratio v — 1/3 [TT| 115]. If the (3D) 
Young's modulus of the sheet is E and its thickness is h, then C s and C& are proportional 
to the 2D Young's modulus and the bending modulus of the sheet: 

C, = ^Eh; C 3 Eh* 

The substrate and actuated layers are assumed to have uniform thicknesses which are 
both equal to h, and the same stretching and bending stiffness constants, for simplicity. In 
the initial configuration, the lower surface of the substrate layer lies in the plane z = —h/2 
and the upper surface lies in the plane z = h/2. The lower surface of the actuated layer 
is bonded to the upper surface of the substrate layer at z = h/2, and the upper surface 



of the actuated layer lies in the z = 3h/2 plane. Throughout the bilayer in its initial, 
flat configuration, there are planar and parallel surfaces of material, lying in the planes of 
constant z between —h/2 and 3h/2. As the bilayer bends, these surfaces of material are no 
longer planar, but are assumed to remain parallel. This assumption underlies the classical 
theories of the bending of thin plates and bilayers (including the continuum versions of ([I]) 
and (|2])), and allows the configuration of the bilayer to be expressed entirely in terms of the 
configuration of the central surface of the substrate (i.e. the material initially lying in the 
plane z = 0). The central surfaces of the substrate and actuated layers may be assumed 
to remain parallel under large deformations, as long as the thickness of the bilayer is small 
compared to its radii of curvature at all points, in all directions. 

The aforementioned triangular mesh represents the central surface of the substrate (i.e. 
the material initially lying in the plane z = 0). The central surface of the actuated layer is 
parallel to that of the substrate, and displaced from it by a distance h in the direction of 
the substrate surface normal n. 

The actuated layer is actuated by setting its equilibrium stretching strain (uniformly here) 
to e a , while the equilibrium strain remains zero for the substrate. Bending causes different 
stretching strains to occur in the central surfaces of the substrate and actuated layers, and 
thereby allows each surface to move closer to its own equilibrium stretching strain. If the 
substrate central surface has curvature k in a certain direction at a point, then the strain in 
the central surface of the actuated layer in that direction is that at the corresponding point 
in the substrate plus hn (see appendix [A]). We consider the central surface of the actuated 
layer to be discretized by a triangular mesh of points, which are those in the substrate mesh 
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plus hn. Then the strain in the actuated layer mesh can be written in terms of that in the 
substrate mesh. Hence only the substrate mesh is required to characterize the full bilayer 
geometry. The lengths of edges in the actuated layer mesh are those of the corresponding 
edges in the substrate mesh, plus the stretching induced by bending in the direction of the 
edge. This stretching is a discrete approximation of the stretching strain Hk, times the 
edge length. In appendix [B] we derive a discrete formula for the edge stretching by using a 
locally-quadratic approximation for the bilayer surface centered at midpoint of each edge. 
The formula is: 

Stretching due to bending « - ip k , (4) 

k=l 

where the angles ip k are those between the normals of the four triangles adjacent to the two 
triangles which share the edge whose stretching is given. These angles are shown schemati- 



cally in figure 16 The formula for the stretching energy of the actuated layer is then: 



1 / hy/3 4 V 
E a ,a = ^Cs I |rj - rj \ H — y^^k- d eq - e a d eq ) . (5) 

i,j \ fc=l 



Equation ^ is similar to ([TJ, but includes the additional stretching due to bending in the 
actuated layer Q, and the equilibrium stretching in the actuated layer e a d eq . In the outer 
sum in E Sta , we omit the contribution from edges which are close to the boundary of the 
mesh, for which not all four of the angles {tfik} are defined. An alternative would be to use 
one-sided approximations for the stretching due to bending near the boundary, which may 
give a higher order of accuracy, at the expense of increased complexity of the scheme. To 
fully realize the increased accuracy in arbitrary initial planar geometries, an unstructured 
triangular mesh may be preferred. However, one of the main benefits of our approach is its 
simplicity. 



The bending energy of the actuated layer is close to that of the substrate layer. In 
appendix [A] we show that the curvature K a in a given direction along the central surface 
of the actuated layer is the same as the curvature of the substrate (at the corresponding 
point and same direction), up to a relative error of Hk. Since Hk <C 1 by assumption in our 
large deformation plate theory, we take the curvature of the actuated layer to be that of the 
substrate. The bending energy of the actuated layer is then the same as §2§. 

The total elastic energy of our model bilayer is then 

E T = E S + E SA + 2E h . (6) 
III. RECTANGULAR BILAYERS 

We now test our scheme for perhaps the simplest bilayers — rectangular bilayers — which 
have been studied previously [TJ [THJ, [T7]. We first consider a rectangular bilayer which is 
nearly square — aspect ratio = a/3/2 = 0.866. The stretching stiffness C s is 8 x 10 4 , and the 
bending stiffness C& is 1. By ([3]), we then have that the thickness of each layer h is 0.01. 
The height and width are of order one (\/3/3 and 2/3 respectively), so we have a thin sheet. 
We set e a to 10/i. At equilibrium, hn is the strain in the actuated layer due to bending in 
the substrate, and this should be comparable to e a . Thus we expect k of about 10. 

We use a mesh with equilibrium edge length 1/60, yielding a discretization with 1661 
points (4983 degrees of freedom) and 4820 edges. We use a large-scale quasi- Newton scheme 
to minimize the discrete energy (|6]), the limited- memory Broyden-Fletcher-Goldfarb-Shanno 
(LM-BFGS) method with a cubic line search. We use an exact formula for the gradient 
of Et which can be computed quickly. An exact gradient formula handles large differences 
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FIG. 1: Bending of a rectangular bilayer discretized with 1661 points and 4820 edges. The stretch- 
ing stiffness C s is 8 x 10 4 , and the bending stiffness C\, is 1. The thickness is h = 0.01, the initial 
height (along y) is "v/3/3, the initial width (along x) is 2/3. The actuation strain e a is lOh. a, 
Configuration of rectangular sheet in equilibrium, with bending along the shorter dimension, b, 
2-norm of energy gradient versus iteration number in LM-BFGS method, leading up to stagnation, 
when the equilibrium is reached, c, Color map of maximum principal curvature at the midpoints of 
edges. Gray edges are too close to the boundary to obtain estimates, d, Direction field showing the 
directions of minimum principal curvature at the midpoints of edges, e, Color map of stretching 
strain in each edge of the substrate central surface. 
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in elastic constants better than a finite difference approximation. The initial guess for the 
minimization routine is the flat state. In figure [T^, we show the equilibrium reached by 
the minimization scheme after it stagnates at about 2500 iterations. The sheet has curled 
along its shorter dimension, yielding a "cigar" shape, one of two equilibria previously found 
for rectangular bilayers [17]. The decrease of the 2-norm of VE? with iteration number is 
shown in panel b, and at stagnation, no further decrease of Et is obtained along the steepest 
descent direction — V-&r- At stagnation, the 2-norm of VE T is 7 x 10~ 3 , and the oo-norm 
is 4 x 10~ 4 . 

Panel c gives a color map showing the larger of the two principal curvatures at the mid- 
points of the edges in the equilibrium configuration in 'a'. Here the bilayer is shown in its 
initial flat state for clarity. The colors of the edges are the same in 'a' and 'c', so the corre- 
spondence between points on the initial and final shape can be seen. We estimate the prin- 
cipal curvatures using the method described in appendix [Cj The estimate can be obtained 
for all edges except those adjacent to the boundary. Over most of the bilayer the maximum 
principal curvature lies between 8.5 and 10.5. Near the boundary, there is more variation. 
Panel d gives a field of line segments showing the directions of minimum-magnitude prin- 
cipal curvature at the edge midpoints. For a developable surface, the minimum-magnitude 
principal curvature is zero at each point. The zero-curvature direction field has integral 
curves that are straight lines, called generators, which cover the surface [H]. The shape 
in panel a (and the other equilbrium shapes in this work) are approximately developable, 
so the lines in panel d approximate the generators of a developable surface. The lines are 
mostly horizontal, showing a cylindrical shape, except near the corners, which curl inward. 
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For cylindrical bending, we can compare our results with an approximate analytical 
continuum solution. Let us assume that the bilayer bends with uniform curvature k in the y 
direction, and zero curvature in x. Let us also assume uniform strains e x and e y in the x and 
y directions. Then, by [TT1 [15] the elastic energy per unit area over the bilayer is uniform 
and equal to 

e T = ^j-C s (3e 2 + 3e 2 + 2e x e y 

AT 

+3(e x - e a ) 2 + 3{e y + hn- e a ) 2 + 2(e x - e a )(e y + hn - e a )) + —C b n 2 . (7) 

In ([7]), the first three of the six terms in the sum multiplying C s correspond to the substrate 
stretching energy, and the second three correspond to the actuated layer stretching energy. 
The last term corresponds to the combined bending energy. The equilibrium, found by 
minimizing ey with respect to e x , e y , and k, is 

e x = y! e y = 0; k = j. (8) 

For the parameters in this simulation, 

e x = 0.05; e y = 0; « = 10. (9) 

In figure [TJ^, most of the curvature values (at edges away from the boundaries) lie between 
8.5 and 10.5, close to the value of k in Q. The curvature is somewhat nonuniform in the 
simulation, unlike in the analytical approximation. Panel e shows the stretching strain in 
each of the edges. For the horizontal edges away from the sheet boundaries, the strain lies 
between 0.045 and 0.055, with an average value close to e x in ([9]). For the other edges (at 
±7r/3 radians from horizontal), the strain lies in the range 0.012-0.016, while that predicted 
by the analytical solution is e x cos 2 n/3 + e y sin 2 n/3 = 0.0125. 
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FIG. 2: Rectangular bilayer bending with the same parameters as in figure [T] and three values 
of equilibrium edge length: 1/30 (a,b), 1/45 (c,d) and 1/60 (e,f). Panels a, c, and e show the 
maximum principal curvatures at the edge midpoints, and panels b, d, and f show the stretching 
strains in each edge. 

In figure [2j we show the results from figure and e together with solutions on two coarser 
meshes, with equilibrium edge lengths of 1/30 (a,b), 1/45 (c,d), and 1/60 (e,f). Panels a, 
c, and e compare maximum principal curvatures using the same color scale. Panels b, d, 
and f compare stretching strains (with nearly the same color scales). The distributions and 
overall magnitudes of these quantities are quite similar from the finest to the coarsest mesh. 
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FIG. 3: Rectangular bilayer bending with the same parameters as in figure [T] and three initial 
aspect ratios: (a) Aspect ratio = 1.15 (initial width = 2/3, initial height = \/3/3 = 0.577), (b) 
Aspect ratio = 2.31 (initial width = 0.8, initial height = v3/5 = 0.346) (c) Aspect ratio = 3.46 
(initial width = 6/7, initial height = V3/7 = 0.247). 

In figure |3j we compare results across three aspect ratios of the initial rectangular shape. 
All three shapes have maximum curvature along the short direction, although for the largest 
aspect ratio (panel c) this is difficult to discern since the strip is much shorter in this 
direction, so it does not show as much overall difference in the z coordinate along this 
direction. Also, it has larger curvature in the long direction than the shapes in panels a and 
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FIG. 4: Bending of an equilateral triangular bilayer with the same elastic parameters as in figure 
[T] The top edge has length 7/9. The values on the color bar indicate the maximum of the 
principal curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color 
map of maximum principal curvatures at each edge midpoint, c, Direction field showing direction 
of minimum principal curvature at the midpoints of edges. 

b. However, for this strip the average curvature along the short direction is 10 times that 
in the longer direction. The average of the curvature in the short direction decreases from 
about 9 for the most square shape (a) to about 8 for the most elongated shape (c). 

IV. GENERAL BILAYERS 

We now consider the bending of a set of bilayers with more general planar shapes. Most 
of the shapes we consider are connected clusters of a small number of equilateral triangles, 
also called "polyiamonds" , a class of polyforms [19] . Such shapes can be represented with 
equilateral triangular meshes of various levels of refinement, without jagged edges (as occur 
for the rectangles in figure |3l for example). 
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FIG. 5: Bending of a V-shaped bilayer with the same elastic parameters as in figure [TJ The bottom 
edge has length 7/9. The values on the color bar indicate the maximum of the principal curvatures 
at each edge midpoint, a, 3D configuration with coloring given in b, Color map of maximum 
principal curvatures at each edge midpoint, c, Direction field showing the directions of minimum 
principal curvature at the midpoints of edges. 

In figure [4] we show the bending of an equilaterial triangular bilayer with the same elastic 
parameters as the rectangular bilayer just considered. The maximum curvature, which 
occurs along the edges, is comparable in magnitude to the maximum curvatures of the 
rectangles in figure |3j but unlike for the rectangles, here there is a large flattened region at 
the center. The minimum-curvature direction field (panel c) shows the triangular symmetry 
of the shape. This equilibrum is not the shape with minimum elastic energy (or global 
equilibrium), however. Cylindrical bending of the triangle is also an equilbrium, and has 
smaller total elastic energy. 

Figure [5] shows the bending of a V-shaped bilayer. Here the two arms of the V bend mainly 
along the longer dimension, in contrast to the rectangular shapes in figure |3j which bend 
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FIG. 6: Bending of a bilayer in the shape of a rhombus with the same elastic parameters as in 
figure [T] The horizontal width is 1/3. The values on the color bar indicate the maximum of the 
principal curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color map 
of maximum principal curvatures at each edge midpoint, c, Direction field showing the directions 
of minimum principal curvature at the midpoints of edges. 

along the shorter dimension. Bending along the longer dimension is also an equilibrium for 
a rectangle, with slightly lower elastic energy than for bending along the shorter dimension 
[T7j . For the V shape, there is a sharp transition to a region of higher curvature near the 
boundary, with bending in the opposite direction. The minimum-curvature direction field 
in panel c shows that near the boundary, there is sharp transition in bending direction. 
Bending is nearly transverse to the boundary near the edges, and parallel to the boundary 
further inside. 

Figure [6] shows the bending of a rhombus. Here the bending is mainly along the shorter 
dimension, akin to the rectangles of figure |3j The direction of bending does not change 
except at the farther pair of opposing tips. Near the edges there are regions of increased 
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FIG. 7: Bending of a tetriamond bilayer with the same elastic parameters as in figure [T] The 
horizontal width is 1/3. The values on the color bar indicate the maximum of the principal 
curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color map of 
maximum principal curvatures at each edge midpoint, c, Direction field showing the directions of 
minimum principal curvature at the midpoints of edges. 

curvature, larger than those in figure |3j 

An oblique strip, a cluster of four equilateral triangles (a tetriamond), is shown in figure 
[7| The bending is essentially cylindrical, and this time along the longer dimension. The 
minimum-curvature direction field is almost orthogonal to the longer edges, but there is a 
noticeable deviation from orthogonality. The curvature is nearly uniform over the middle 
part of the bilayer, and varies more near the tips with acute angles. The final shape maintains 
the 180-degree rotational symmetry of the initial shape. 

A different tetriamond is shown in figure [8] The bending is again essentially cylindrical, 
but this time along the shorter dimension. Now the minimum-curvature direction field is 
essentially orthogonal to the upper and lower edges, giving a shape with bilateral symmetry. 
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FIG. 8: Bending of a tetriamond bilayer with the same elastic parameters as in figure [T] The 
horizontal width is 1/3. The values on the color bar indicate the maximum of the principal 
curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color map of 
maximum principal curvatures at each edge midpoint, c, Direction field showing the directions of 
minimum principal curvature at the midpoints of edges. 

A pentiamond, or hexagon with a wedge removed, is shown in figure [9j The 3D config- 
uration has a bilateral symmetry inherited from the 2D shape. The bending is not quite 
cylindrical. The minimum-curvature direction field gradually rotates, moving around the 
central corner. At the central corner, the curvature is significantly increased, perhaps indi- 
cating a smoothed version of a conical singularity. 



Figure 10 shows the bending of a different pentiamond. Again, the curvature is largest 



at the concave angle. The minimum-curvature directions are mainly vertical on the right 
side of panel c, and transition to an oblique angle on the left side, yielding an approximate 
partition into two regions of bending. Overall, the bending is mainly along the longest 



dimension of the shape. 
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FIG. 9: Bending of a pentiamond bilayer with the same elastic parameters as in figure [TJ The 
horizontal width is 1/3. The values on the color bar indicate the maximum of the principal 
curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color map of 
maximum principal curvatures at each edge midpoint, c, Direction field showing the directions of 
minimum principal curvature at the midpoints of edges. 



A strip-like pentiamond is shown in figure 11 similar to the tetriamond of figure [7j Unlike 



that shape, here the initial and final shapes have bilateral symmetry. Similarly to the shape 
of figure [7j the bending is mainly along the longer dimension, and curvature varies most 
near the two ends. 



In figure 12, another pentiamond bilayer is shown. The shape has no symmetries, and 
the overall bending is strip-like, with the minimum-curvature directions mainly orthogonal 
to the longest linear dimension of the object. 



The last of the four pentiamonds is shown in figure 13 , another shape without symmetries 



The bending is most similar to that of figure 10 with the minimum-curvature direction field 



again split mainly into two regions, with some amount of convergence at the concave angle. 
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FIG. 10: Bending of a pentiamond bilayer with the same elastic parameters as in figure[TJ The top 
edge has length 1/2. The values on the color bar indicate the maximum of the principal curvatures 
at each edge midpoint, a, 3D configuration with coloring given in b, Color map of maximum 
principal curvatures at each edge midpoint, c, Direction field showing the directions of minimum 
principal curvature at the midpoints of edges. 

The bending is mainly along the longest dimension. 



We consider only two of the twelve hexiamonds. In figure 14 , a hexagonal bilayer is shown 



The bending is roughly cylindrical, and nearly along a line connecting two opposite vertices. 
The curvature distribution has a bilateral symmetry, with the strongest intensification of 
curvature at the sides parallel to the bending direction. 



A second hexiamond, with rotational and reflectional symmetries, is shown in figure 15 



The bending is mainly cylindrical here also, but there is a deviation in bending direction 
near the 60-degree corners, similar to that which has occured in most but not all of the 
60-degree corners in the bilayers previously discussed. These corners curl inward toward the 



center of the shape. 
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FIG. 11: Bending of a pentiamond bilayer with the same elastic parameters as in figure [I] The 
bottom edge has length 1/4. The values on the color bar indicate the maximum of the principal 
curvatures at each edge midpoint, a, 3D configuration with coloring given in b, Color map of 
maximum principal curvatures at each edge midpoint, c, Direction field showing the directions of 
minimum principal curvature at the midpoints of edges. 

We have presented just a single equilibrium shape for a variety of bilayer shapes, most 
corresponding to simple polyiamonds. These shapes show how the simple cylindrical bending 
of figures [T] and [3] are modified for oblique shapes and shapes with different or no symmetries. 
These examples may be useful test cases from which to abstract the general dependence of 
bilayer bending on initial shape. 

V. CONCLUSION 

We have presented a simple discrete formula for the elastic energy of a bilayer, in terms 
of the geometry of a single mesh representing the central surface of the substrate layer. The 
formula allows for fast simulations of bending bilayers with general geometries. We have 
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FIG. 12: Bending of a pentiamond bilayer with the same elastic parameters as in figure [TJ The top 
edge has length 2/3. The values on the color bar indicate the maximum of the principal curvatures 
at each edge midpoint, a, 3D configuration with coloring given in b, Color map of maximum 
principal curvatures at each edge midpoint, c, Direction field showing the directions of minimum 
principal curvature at the midpoints of edges. 

found good agreement between our simulations of a rectangular bilayer, and an approximate 
analytical solution. 

Simulations of a set of shapes composed of clusters of equilateral triangles (polyiamonds) 
present some generic behaviors of bilayers which may lead to an understanding of the rela- 
tionship between bending directions and bilayer shapes. Some of the bilayers show cylindrical 
bending in the same direction over the entire bilayer, and often along either the shortest 
or longest linear extent of the initial shape. Slight deviations in bending direction occur 
sometimes but not always near acute-angled corners of the bilayers, where the shape bends 
inwards towards the center. Larger variations in curvature occur near edges and corners, and 



the curvature is particularly intensified near obtuse-angled corners, as in the pentiamond 
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FIG. 13: Bending of a pentiamond bilayer with the same elastic parameters as in figure[TJ The top 
edge has length 1/3. The values on the color bar indicate the maximum of the principal curvatures 
at each edge midpoint, a, 3D configuration with coloring given in b, Color map of maximum 
principal curvatures at each edge midpoint, c, Direction field showing the directions of minimum 
principal curvature at the midpoints of edges. 

which is a hexagon with a sector removed. In this shape and some others, the direction 
of bending changes gradually, leading to bending which is somewhat conical rather than 
cylindrical. Some other shapes are nearly partitioned into two regions, with bending along 
distinct directions in each region. 

Because the bilayers are thin, the equilibria are all close to developable surfaces. Unlike 
thin homogeneous sheets, in bilayers there is always considerable in-plane stretching energy, 
in the direction of small curvature. The equilibria we have found preserve many of the initial 
symmetries of the flat bilayer. Performing our simulations with other initial guesses for the 
equilibria can produce different and possibly less symmetrical equilibria. 
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FIG. 14: Bending of a hexagonal bilayer with the same elastic parameters as in figure [TJ The sides 
have length 1/3. The values on the color bar indicate the maximum of the principal curvatures at 
each edge midpoint, a, 3D configuration with coloring given in b, Color map of maximum principal 
curvatures at each edge midpoint, c, Direction field showing the directions of minimum principal 
curvature at the midpoints of edges. 

Appendix A: Extension due to bending 



Let us consider a point X on the central surface of the substrate. Let the curvature at 
Xo in a direction v be k, and the radius of curvature be R = 1/k. Consider the intersection 
of the surface with the plane spanned by v and the normal to the surface, n. This is a 
curved line passing through Xo. A short segment of length I near Xo is approximated by an 
arc of a circle tangent to the surface at X with radius R. The set of corresponding points 
on the central surface of the actuated layer is a parallel curve, approximated by a concentric 
arc with radius R + h, and length / + 61. Because the arcs are concentric, 

^ = *±* and thus S 4 = h K , (Al) 
/ R I 



23 



FIG. 15: Bending of a hexiamond with the same elastic parameters as in figure [TJ The four shorter 
sides have length 1/3 and the two longer sides have length 2/3. The values on the color bar indicate 
the maximum of the principal curvatures at each edge midpoint, a, 3D configuration with coloring 
given in b, Color map of maximum principal curvatures at each edge midpoint, c, Direction field 
showing the directions of minimum principal curvature at the midpoints of edges. 

so the strain in the actuated layer is hn plus that in the substrate. Furthermore, the ratio 
of the curvature of the actuated layer central surface K a to that of the substrate k is 

= 5TS = TTfwi = 1 - hK + ° (i hK ) 2 ) ■ (A2) 

Appendix B: Discrete formula for extension due to bending 

In order to derive equation Q, we consider the discretization in the vicinity of a point 
on the central surface of the substrate, which is assumed to be a smooth surface. Without 
loss of generality, we can assume that the point is at the origin in (x, y, z), and that, up to 
quadratic order, the central surface of the subtrate is given by 

z ( x , y) = \.Kx 2 + \k y y 2 - ( B1 ) 
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FIG. 16: Schematic diagram of a set of triangles showing the notation used in deriving the formula 
for the strain in edge X1-X2 due to bending. Here the amount of bending is exaggerated so the 
differences between vectors are visible. 

An arbitrary smooth surface can be brought into this form (up to quadratic order), in the 
vicinity of a point on the surface, by translating the point to the origin and rotating the 
surface about the point so that the directions of principal curvature lie along the x and y 
axes. 

We now assume that the origin coincides with the midpoint of an edge in our triangulation 
of the substrate central surface. We determine the stretching of this edge in the actuated 



layer relative to that in the substrate, under the bending of the substrate given by (Bl). 



In figure 16 we show a portion of the triangular mesh consisting of the two triangles which 



bound the X1-X2 edge, and the four other triangles which share an edge with these triangles. 
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We can approximate the curvatures of the surface simply in terms of the angles between the 
normals to these triangles. In general the mesh has an arbitrary orientation with respect to 



the directions of principal curvature (here, the x and y axes). Thus in figure 16 we assume 



an arbitrary angle <fi between the x-axis and the projection of the X1-X2 edge on the x-y 



plane. We first assume that the x-y coordinates of the eight points in figure 16 lie on the 



undeformed equilateral triangular lattice, which is a good approximation for small strains: 



Xl,: 



(X 3ia5 ,X 3) „ 



(X4 >x , X4 j?/ 



(X5 iX , Xs^ 



(X6,a;, X6,j, 
(Xg )X , Xg^ 



d, 



cos0, X . , ; =:- — sni < 



2 



— (Xi )X ,X 

(Xl >X , Xl f y 



l,y) 



(Xi )X ,X 



1,2/ 



(Xg,!, Xg^ 



+ d eq ( cos 



7T 



, sin 



(cos 0, sin ( 



— deq (cos (/>, sin ( 



+ d e9 (cOS (4> + 



— c? e9 (cos 0, sin ( 



G? eg (cos 0, sin ( 



7T 



37' 



sin 



7T 



7T 



(B2) 
(B3) 
(B4) 
(B5) 
(B6) 
(B7) 
(B8) 
(B9) 



The z coordinates of Xi-Xg are found by inserting the x-y coordinates in (B2)-(B9) into 



(Bl ). The points which correspond to Xi-X 2 in the central surface of the actuated layer are 



Xi = Xi + hni 
X 2 = X 2 + hn 2 



(BIO) 
(Bll) 
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where ni and ri2 are the unit normal vectors to the surface (Bl ) at Xi and X2, respectively: 



ni = f-k x ^ cos0, -ky^ sin 0> !^ (l + O {dl q )) 
n 2 = \ cos< ^' k v~Y sin< ^' X ) ( X + ( d eq)) ■ 



(B12) 
(B13) 



We may now determine the length of the edge Xi-X 2 relative to d eq : 



Xx — Xj 



4 ? = -hd eq {k x cos 2 + fcj, sin 2 <f>)(l + (d 2 eq ) ) (1 + O (h)) . (B14) 



If the corresponding edge Xx-X 2 in the substrate has a length d different from d eq , then 



(B14) still holds but with d in place of d eq . Thus it gives the difference in stretching between 



the actuated layer central surface and the substrate at corresponding edges. 



We now show that (B14) is approximately — /ia/3 times the average of the angles between 
the normals to the four triangles adjacent to the two triangles which share the X1-X2 edge. 
We define vectors normal to triangles 1-6 



n x = 


(X! 


-x 3 ) 


X 


(X 4 


-x 3 ) 


(B15) 


n 2 = 


(X 4 


-x 2 ) 


X 


(X 2 


-Xx) 


(B16) 


n 3 = 


(x 4 


-x 5 ) 


X 


(X 2 


-x 5 ) 


(B17) 


n 4 = 


(X 7 


-x 6 ) 


X 


(X! 


-x 6 ) 


(B18) 


n 5 = 


(x 7 


-xo 


X 


(x 2 


-Xx) 


(B19) 


n 6 = 


(X 7 


-x 2 ) 


X 


(X 8 


-x 2 ) 


(B20) 



The angle between n 1 and n 2 is 



(pi = Arccos j, — - 

\ Til 



n 1 ■ n 2 



(B21) 
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Using ( B15|B16 ) we define a function a\ by 

^n 1 ■ n 2 = d%{\ + d 2 eqai (k x , k y , 0)), (B22) 



a\ has terms which are quadratic in k x and k y and include trigonometric functions of <j). We 



also define b\ and C\ for terms in the denominator of (B21): 

4 



3 
4 
3 



n 1 f = d i eq {l + d 2 eq b l {k x ,k y ,<f))) (B23) 



L n 2 f = ^(1 + ^(^,^,0)) (B24) 



Thus 



In 1 !! Iln 2 | 



nl ' n2 l + + (B25) 



We define a 2 and 62 in terms of n 2 and n 

4 



3 
4 
3 



n 3 ■ n 2 = d 4 eq (l + d 2 eq a 2 (k x , ky, 0)) (B26) 



|n 3 || 2 = <(l + d 2 eq b 2 (k x ,k y ,^)) (B27) 



n 2 • n 3 



The angle between n 2 and n 3 is 

to = A ^cos ^1^,111^11 j (B28) 
We expand Arccos in a Taylor series 

Arccos(l -x) = V2x + 0(x 3/2 ). (B29) 

Thus 



\{to + to) = ^d eg f \j- ai + * + * + x /- a2 + ] (l + O «)) (B30) 
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Using trigonometric identities one can show that 



61 + Ci / 6 2 + Ci 



f fca, cos 2 4> + k y sin 2 



(B31) 



Therefore by (B14), 



- hV3 



<Pi + ¥2 



(l + 0«)) (1 + O0O) 



d, 



eq- 



(B32) 



A more symmetrical formula is obtained by using the angles 



(fz = Arccos 



4 5 
n • n 



|n 4 || ||n 5 | 



if 4 = Arccos 



n 5 • n 6 



|n 5 ||||n 6 | 



(B33) 



By symmetries of the lattice (B2)-(B9) and the quadratic surface, ^3 = ^2 and y? 4 = <p\ 



Thus a more symmetrical alternative to (B32) is 



- hV3 



4 



(i + o«)) (1 + 00) 



Xi — X- 



-d, 



eq- 



(B34) 



(B34) is more accurate than (B32) for surface approximations which are of higher than 



quadratic order, in which case the principal curvatures have a nonzero gradient at the origin. 
Let us now assume that the substrate central surface is strained from the undeformed 



triangular mesh, given by (B2)-(B9), with a nonzero in-plane strain tensor c^,-. The deriva- 



tion leading to (B34) can be repeated, but now there is an additional factor of 1 + 0(||a;| 



multiplying the left side of (B34). This error is of the same order as errors already present 



in the single-plate model represented by ([!]) and (J2J) alone. These errors are due to the use 
of only the linear part of the strain tensor in infmitessimal strain theory, and also due to 
the discrete approximation ([!]) (see [TT]). 
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Appendix C: Estimate of principal curvatures at a point 



We now use the framework from appendix[B]to obtain a discrete estimate for the principal 
curvatures at the midpoint of an edge in our mesh. We consider the edge shown in figure 



16, oriented at an angle (p with respect to the x axis, and ir/2 — <p with respect to the y axis. 



The x and y axes are again assumed to be the directions of principal curvatures, and near 



the edge midpoint the surface is (Bl ). We can use (B34) to estimate the curvature k\ in the 



direction of the edge X]-X 2 (i.e. along the angle in figure 16): 



hk\d, 



eq 



<fi + ¥2 + V?3 + (Pa 
4 



(CI) 



We can similarly estimate the curvatures along the directions (f> ± ir/3. Let fc 2 denote the 



curvature in the direction + 7r/3, at the midpoint of X!-X 2 . We apply formula (CI ) to the 



edges X2-X7 and X1-X4 which have angle <fi + tt/3, and take the average as an estimate of 



the curvature at the midpoint of Xi-X 2 . We take the average of (CI) applied to the edges 



Xx-X 7 and X 2 -X 4 to estimate k 3 , the curvature in direction <fi — ir/3. 

The curvatures k\-k 3 are related to the principal curvatures k x and k y by 



k\ = k x cos 2 <p + k y sin 2 0. 

fc 2 = k x cos 2 \ $ Jr ~^j ^ k y sin 2 

n \ 7-2 

-)+k y sin 



+ 



k 3 = k x cos 2 ^ 



(C2) 
(C3) 
(C4) 



Given our estimates of ki-k 3 at the midpoint of a given edge, we regard (C2)-(C4) as a system 



of three nonlinear equations to be solved for three unknowns: the principal curvatures k a 



and k v , and the angle <fi which orients the mesh with respect to the directions of principal 
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curvature. The equations are simpler in terms of the alternate set of three variables 



{k y , u, A}, u = cos 2 4>, A = k x — k y . 



(C5) 



The solutions are 



u 



1 (2fci - fc 2 - fcg) sign(A;2 - fcg) 

2 2 V /3(A;3-A; 2 ) 2 + (2A; 1 -A; 2 -A;3) 2: 



.4 



VM1 



u 



k y = k 1 - Au. (C6) 



Using (C6) and (C2)-(C4) we can estimate the principal curvatures at the midpoints of all 



the edges in our mesh except those near the boundary, for which not all the edges needed 



to estimate k±, k 2 , and k 3 exist. 
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